Quantum dot self consistent electronic structure and the Coulomb blockade 
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We employ density functional theory to calculate the self consistent electronic structure, free 
energy and linear source-drain conductance of a lateral semiconductor quantum dot patterned via 
surface gates on the 2DEG formed at the interface of a GaAs — AlGaAs heterostructure. The 
Schrodinger equation is reduced from 3D to multi-component 2D and solved via an eigenfunction 
expansion in the dot. This permits the solution of the electronic structure for dot electron number 
N ~ 100. We present details of our derivation of the total dot-lead-gates interacting free energy 
in terms of the electronic structure results, which is free of capacitance parameters. Statistical 
properties of the dot level spacings and connection coefficients to the leads are computed in the 
presence of varying degrees of order in the donor layer. Based on the self-consistently computed 
free energy as a function of gate voltages, Vi, and N, we modify the semi-classical expression for 
the tunneling conductance as a function of gate voltage through the dot in the linear source-drain, 
Coulomb blockade regime. Among the many results presented, we demonstrate the existence of a 
shell structure in the dot levels which (a) results in envelope modulation of Coulomb oscillation peak 
heights, (b) which influences the dot capacitances and should be observable in terms of variations 
in the activation energy for conductance in a Coulomb oscillation minimum, and (c) which possibly 
contributes to departure of recent experimental results from the predictions of random matrix theory. 

PACS numbers: 73.20.Dx,73.40.Gk,73.50. Jt 



I. INTRODUCTION 

Study of the Coulomb blockade and charging effects in 
the transport properties of semiconductor systems is pe- 
culiarly suitable to investigation through self-consistent 
electronic structure techniques. While the orthodox 
theorya, in parameterizing the energy of the system in 
terms of capacitances, is strongly applicable to metal 
systems, the much larger ratio of Fermi wavelength to 
system size, Xf/L, in mesoscopic semiconductor devices, 
requires investigation of the interplay of quantum me- 
chanics and charging. 

In the first step beyond the orthodox theory, the "con- 
stant interaction" model of the Coulomb blockade supple- 
mented the capacitance parameters, which were retained 
to characterize the gross electrostatic contributions to the 
energy, with non-interacting quantum levels of the dots 
and leads of the mesoscopic deviceatl. This theory was 
successful in explaining some of the fundamental features, 
specifically the periodicity, of Coulomb oscillations in the 
conductance of a source-dot-drain-gate system with vary- 
ing gate voltage. Other effects, however, such as varia- 
tions in oscillation amplitudes, were not explained. 

In this paper we employ density functional (DF) the- 
ory to compute the self-consistently changing effective 
single particle levels of a lateral GaAs — AlGaAs quan- 
tum dot, as a function of gate voltages, temperature T, 
and dot electron number Na. We also compute the total 
system free energy from the results of the self-consistent 
calculation. We are then able to calculate the device con- 



ductance in the linear bias regime without any adjustable 

parameters. Here we consider only weak (< 0.1 T) mag- 
netic fields in order to study the effects of breaking time- 
reversal symmetry. We will present results for the edge 
state regime in a subsequent publicationO. 

We include donor layer disorder in the calculation and 
present results for the statistics of level spacings and par- 
tial level widths due to tunneling to the leads. Recently 
we have employed Monte-Carlo variable range hopping 
simulations to consider the effect of Coulomb regulated 
ordering of ions in the donor layer^n the mode character- 
istics of split-gate quantum wireM. The results of those 
simulations are here applied to quantum dot electronic 
structure. 

A major innovation in this calculation is our method 
for determining the two dimensional electron gas (2DEG) 
charge density. At each iteration of the self-consistent 
calculation, at each point in the x — y plane we determine 
the subbands e n (x,y) and wave functions £n V ( z ) m the 
z (growth) direction. The full three dimensional density 
is then determined by a solution of the multi-component 
2D Schrodinger equation and/or 2D Thomas- Fermi ap- 
proximation. 

Among the many approximation in the calculation are 
the following. We use the local density approximation 
(LDA) for exchange-correlation (XC), specifically the pa- 
rameterized form of Stern and Das SarmaQ. While the 
LDA is difficult to justify in small (N ~ 50 - 100) 
quantum dots it is empirically known to give good re- 
sults in atomic and molecular systems where the density 
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is also changing appreciably on the scale of the Fermi 
wavelengtba. 

In reducing the 3D Schrodinger equation to a multi- 
component 2D equation we cutoff the expansion in sub- 
bands, often taking only the lowest subband into account. 
We also cutoff the wavefunctions by placing another ar- 
tificial AlGaAs interface at a certain depth (typically 
o 

200 A) away from the first interface, thereby ensuring 
the existence of subbands at all points in the x — y plane. 
Generally the subband energy of this bare square well is 
much smaller than the triangular binding to the interface 
in all but those regions which are very nearly depleted. 

The dot electron states in the zero magnetic field 
regime are simply treated as spin degenerate. For B ^ 
an unrenormalized Lande g-factor of —0.44 is used. We 
employ the effective mass approximation uncritically and 
ignore the effective mass difference between GaAs and 
AlGaAs (m* — 0.067 m ). Similarly we take the back- 
ground dielectric constant to be that of pure GaAs (k — 
12.5) thereby ignoring image effects (in the AlGaAs). 
We ignore interface grading and treat the interface as a 
sharp potential step. These effects have been treated in 
other calculations of self-consistent electronic structure 
for GaAs — AlGaAs devicesQ and have generally been 
found to be small. 

We mostly employ effective atomic units wherein 
I Ry* = m*e 4 /2h 2 K 2 w 5.8 meV and f a% = 

H 2 n/m*e 2 w 100 A- 

The structure of the paper is as follows. In section II 
we first dismiss the calculation of the electronic structure, 
focusing on those features which are new to our method. 
Further subsections then consider the treatment of dis- 
crete ion charge and disorder, calculation of the total dot 
free energy from the self-consistent electronic structure 
results, calculation of the source-dot-drain conductance 
in the linear regime and calculation of the dot capac- 
itance matrix. Section III provides new results which 
are further subdivided into basic electrostatic properties, 
properties of the effective single electron spectra, statis- 
tics of level spacings and widths and conductance in the 
Coulomb oscillation regime. Section IV summarizes the 
principal conclusions which we derive from the calcula- 
tions. 



II. CALCULATIONS 

A. Quantum dot self-consistent electronic structure 

We consider a lateral quantum dot patterned on a 
2DEG heterojunction via metallic surface gates (Fig. [I]). 
At a semiclassical level, other gate geometries, such as 
a simple point contact or a multiple dot system, can be 
treated with the same methodQ'El. However, a full 3D 
solution of Schrodinger's equation, even employing our 
subband 
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FIG. 1. Schematic of device used in calculation. The 
z-subband structure throughout the plane are calculated at 
each iteration of the self-consistency loop. Most results pre- 
sented with gate variation assume that both the upper and 
lower pins of the relevant gate are simultaneously varied. 

expansion procedure for the z direction, is only 
tractable in the current method when a region with a 
small number of electrons (N < 100) is quantum me- 
chanically isolated, such as in a quantum dot. 



1. Poisson equation and Newton's method 

In principal, a self-consistent solution is obtained by 
iterating the solution of Poisson's equation and some 
method for calculating the charge density (see following 
sections II. A. 2 and II. A. 3). In practise, we follow Kumar 
et aild and use an A/"-dimensional Newton's method for 
finding the zeroes of the functional F((f>) = A.-(f>+p((f>)+q; 
where the potential, (pi, and density pi, on the Af discrete 
lattice sites (Af ~ 100, 000) are written as vectors, <f> and 
p. The vector q represents the inhomogeneous contri- 
bution from any Dirichlet boundary conditions, A is the 
Laplacian (note that here it is a matrix, not a differential 
operator), modified for boundary conditions. Innovations 
for treating the Jacobian dpi/d(j>j beyond 3D Thomas- 
Fermi, and for rapidly evaluating the mixing parameter t 
(see RefJiB) are discussed below. 

The Poisson grid spans a rectangular solid and hence 
the boundary conditions on six surfaces must be sup- 
plied. Wide regions of the source and drain must be 
included in order to apply Neumann boundary condi- 
tions on these (x = constant) interfaces, so a non-uniform 
mesh is essential. It is also possible to apply Dirichlet 
boundary conditions on these interfaces using the un- 
gated wafer (one dimensional) potential profile calculated 
off-linetil. In this case, failure to include sufficiently wide 
lead regions shows up as induced charge on these surfaces 
(non- vanishing electric field). To keep the total induced 
charge on all surfaces below 0.5 electron, lead regions of 
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~ 5 \im are necessary, assuming a surface gate to 2DEG 

o 

distance (i.e. AlGaAs thickness) of 1000 A- In other 
words we need an aspect ratio of 50 : 1. We note that 
we ignore background compensation and merely assume 
that the Fermi level is pinned at some fixed depth ("zoo" 
^2.5 urn) into the GaAs at the donor level. The donor 
energy for GaAs is taken as 1 Ry* below the conduction 
band. In the source and drain regions, the potential of 
the 2DEG Fermi surface is fixed by the desired (input) 
lead voltage. 

We apply Neumann boundary conditions at the y = 
constant surfaces. The z = surface of the device has 
Dirichlet conditions on the gated regions (voltage equal 
to the relevant desired gate voltage) and Neumann con- 
ditions, dcf>/dn = 0, elsewhere. This is equivalent to the 
"frozen surface" approximation ofll3, further assuming a 
high dielectric constant for the semiconductor relative to 
air. Further discussion of this semiconductor-air bound- 
ary condition can be found in Ref.E3. 



2. Charge density, quasi-2D treatment 

The charge density within the Poisson grid (i.e. not 
surface gate charge) includes the 2DEG electrons and 
the ions in the donor layer. The treatment of discrete- 
ness, order and disorder in. the donor ionic charge pi on 
has been discussed in Ref.cl in regards to quantum wire 
electronic structure. Some further relevant remarks are 
made below in section II. B. 

As noted above, we take advantage of the quasi-2D na- 
ture of the electrons at the GaAs — AlGaAs interface to 
simplify the calculation for their contribution to the total 
charge. Given <j), we begin by solving Schrodinger's Eq. 
in the z-direction at every point in the x — y plane, 

h + V B {z) + e<j>(x,y, z)]C v (z) = e n (x, y)Q y (z) (1) 

where Vb(z) is the potential due to the conduction band 
offset between GaAs and Al x Ga\- x As. We generally 
employ fast Fourier transform with 16 or 32 subbands. 

In order that there be a discrete spectrum at each point 
in the x — y plane, it is convenient to take Vb(z) as a 
square well potential (Fig. |l]). That is, we effectively 
cutoff the wave function with a second barrier, typically 


200 A from the primary interface. In undepleted regions 
the potential is still basically triangular and only the tail 
of the wave function is affected. However, near the bor- 
der between depleted and undepleted regions the arti- 
ficial second barrier will introduce some error into the 
electron density. This is because as a depletion region is 
approached, the binding electric field at the 2DEG inter- 
face (slope of the triangular potential) reduces, in addi- 
tion to the interface potential itself rising. Consequently, 
all subbands become degenerate and near the edge elec- 
trons are three dimensional^. We have checked that this 



departure from interface confinement, and in general in- 
plane gradients of £n V (z) contribute negligibly to quan- 
tum dot level energies. However, theoretical descriptions 
of 2DEG edges commonly assume perfect confinement of 
electrons in a plane. In particular the description of edge 
excitations in the quanturp,Hall effect regime in terms 
of a chiral Luttinger liquidEj may be complicated in real 
samples by the emergence of this vanishing energy scale 
and collective modes related to it. 

Assuming only a single z-subband now and dropping 
the index n, we determine the charge distribution in the 
x — y plane from the effective potential e(x, y), employing 
a 2D Thomas-Fermi approximation for the charge in the 
leads and solving a 2D Schrodingcr equation in the dot. 
In order that the dot states be well defined, the QPC sad- 
dle points must be classically inaccessible. (If this is not 
the case it is still possible to use a Thomas- Fermi approx- 
imation throughout the plane for the charge densitytm). 
In the dot, the density is determined from the eigenstates 
by filling states according to a Fermi distribution cither 
to a prescribed "quasi- Fermi energy" of the dot, or to a 
fixed number of electrons. It has been pointed out that a 
Fermi distribution for the level occupancies in the dot is 
an inaccurate approximation to the correct grand canon- 
ical ensemble distribution^. Nonetheless, for small dots 
(A ~ 15) Jovanovic et a/.0 have shown that, regard- 
ing the filling factor, the discrepancy between a Fermi 
function evaluation and that of the full grand canonical 
ensemble is ~ 5% at half filling and significantly smaller 
away from the Fermi surface. As N increases the dis- 
crepancy should become smaller. 



3. Solution of Schrodinger's equation in the dot 

To solve the effective 2D Schrodinger's equation in the 
dot, 

(-V 2 + e (x))/(x) =J B/(x) (2) 

we set the 2D potentials throughout the leads to their val- 
ues at the saddle points, thereby ensuring that the wave 
functions decay uniformly into the leads. Thus the energy 
of the higher lying states will be shifted upward slightly. 
In seeking a basis in which to expand the solution of 
Eq. |^ we must consider the approximate shape of the 
potential. The quantum dots which we model here are 
lithographically approximately square in shape. However 
the potential at the 2DEG level and also the effective 2-D 
potential e(r, 9), (now in polar coordinates) are to lowest 
order azimuthally symmetric. The radial dependence of 
the potential is weakly parabolic across the center. Near 
the perimeter higher order terms become important (cf. 
figure |b and Eq. |H). 

As the choice of a good basis is not completely clear, 
we have tried two different sets of functions: Bessel func- 
tions and the so-called Darwin-Fock (DF) statesE-3. The 
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details of the solution for the eigenfunctions and eigen- 
values differ significantly whether we use the Bessel func- 
tions or the DF states. The Bessel function case is largely 
numerical whereas the DF functions together with poly- 
nomial fitting of the azimuthally symmetric part of the 
radial potential allow a considerable amount of the work 
to be done analytically. Further, neither of the two bases 
comes particularly close to fitting the somewhat eccentric 
shape of the actual dot potential. It is therefore gratify- 
ing that comparing the eigenvalues determined from the 
two bases when reasonable cutoffs are used, we find for 
up to the 50 th eigenenergy agreement to three significant 
figures, or to within roughly 5 micro eV. 



4- Summary and efficiency 

To summarize the calculation, we begin by choosing 
the device dimensions such as the gate pattern, the ion- 
ized donor charge density and its location relative to the 
2DEG, the aluminum concentration for the height of the 
barrier, and the thickness of the AlGaAs layer. We con- 
struct non- uniform grids 'm x, y and z that best fit the 
device within a total of about 10 5 points. Gate volt- 
ages, temperature, source-drain voltages, and either the 
electron number N or the quasi-Fermi energy of the dot 
are inputs. The iteration scheme begins with a guess 
of . The TD Schrodinger equation is solved at each 
point in the x — y plane and an effective 2-D potential 
e(x,y) for one or at most two subbands is thereby de- 
termined. Taking iCn 1 ^)! 2 for the z-dependence of the 
charge density, we compute the 2D dependence in the 
leads using a 2D Thomas-Fermi approximation and in 
the dot by solving Schrodinger's equation and filling the 
computed states according to a Fermi distribution. We 
compute F{(f>^), which is a measure of how far we are 
from self-consistency, and solve for Sip, the potential in- 
crement, using a mixing parameter t. This gives the next 
estimate for the potential (fy- 1 ' . The procedure is iterated 
and convergence is gauged by the norm of F. 

In practise there are many tricks which one uses to 
hasten (or even obtain !) convergence-,.-Eirst, we use a 
scheme developed by Bank and RoseHL3 to search for 
an optimal mixing parameter t. Repeated calculation of 
Schrodinger's equation, which is very costly, is in prin- 
ciple required in the search for t. Far from convergence 
the Thomas-Fermi approximation can be used in the dot 
as well as the leads. Nearer to convergence we find that 
diagonalizing t 5<fi in a basis of about ten states near 
the Fermi surface, treating the charge in the other filled 
states as inert, is highly efficient. Periodically the full 
solution of Schrodinger's equation is employed to update 
the wave functions. 

The wave function information is also used to make 
a better estimate of dpi/d(fii. The 3D Thomas-Fermi 
method for estimating this quantity does not account for 
the fact that the change in density at a given grid point 



will be most strongly influenced by the changes in the 
occupancies of the partially filled states at the Fermi sur- 
face. Thus use of these wave functions greatly improves 
the speed of the calculation. 

B. Disorder 

Evidence of Coulombic ordering of the donor charge 
in a modulation doping layer adjacent to a 2DEG has 
recently accumulated!!! When the fraction T of ionized 
donors among all donors is less than unity, redistribution 
of the ionized sites throushJjiopping can lead to ordering 
of the donor layer chargeEMl 

In this paper we consider the effects of donor charge 
distribution on the statistical properties of quantum dot 
level spectra, in particular the unfolded level spacings, 
and on the connection coefficients to the leads T p of the 
individual states (see below). These dot properties are 
calculated with ensembles of donor charge which range 
from completely random (identical to T = 1, no ion re- 
ordering possible) to highly ordered [T ~ 1/10). For a 
discussion of the glass-like properties of the donor layer 
and the Monte-Carlo variable range hopping calculation 
which is used to generate ordered ion ensembles, see 
Refs.B and 2 ! 

Note that hopping is assumed to take place at tempera- 
tures (~ 160 K) much higher than the sub-liquid Helium 
temperatures at which the dot electronic structure is cal- 
culated. Thus the ionic charge distributions generated in 
the Monte-Carlo calculation are, for the purposes of the 
2DEG electronic structure calculation, considered fixed 
space charges which are specifically not treated as being 
in thermal equilibrium with the 2DEG. 

The region where the donor charge can be taken as dis- 
crete is limited by grid spacing and hence computation 
time. In the wide lead regions and wide region lateral to 
the dot the donor charge is always treated as "jellium." 
Also, to serve as a baseline, we calculate the dot structure 
with jellium across the dot region as well. We introduce 
the term "quiet dot" to denote this case. 

C. Free energy 

To calculate the total interacting free energy we begin 
from the semi-classical expression 

F({n p }, Qi, V t ) = E P n P e° p + \ £f Q l V i 

-Ei^otfdtVimit) (3) 

where n p are the occupancies of non-interacting dot en- 
ergy levels e p ; Qi and Vj are the charges and voltages of 
the M distinct "elements" into which we divide the sys- 
tem: dot, leads and gates. Ii are the currents supplied 
by power supplies to the elements. 

The self-consistent energy levels for the electrons in 
the dot are e p =< ip p | -V 2 + V B (z) + e<j>{r) \ip p >. A 
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sum over these levels double counts the electron-electron 
interaction. Thus, for the terms in Eq. || relating to the 
dot, we make the replacement: 



E p n P £ p + \QdotVdot -> E P n P £ P 
| / drp dot (r)(j)(r) + \ J drp io „(r)</>(r) 



(4) 



where Pdot(?) refers only to the charge in the dot states 
and /?ion(r) refers to all the charge in the donor layer. 

We have, |— - demonstratedEfE 2 ] that previous 
investigationsuS had failed to correctly include the work 
from the power supplies, particularly to the source and 
drain leads, in the energy balance for tunneling between 
leads and dot in the Coulomb blockade regime. Here, we 
assume a low impedance environment which allows us to 
make the replacement: 



(5) 



i^dot 



i^dot 



i^dot 



The charges on the gates are determined from the gra- 
dient of the potential at the various surface regions, the 
voltages being given. Including only the classical electro- 
static energy of the leads, the total free energy isQ: 

F({n p }, N, V) = J2 P ripSp - § / drp dot {r)(f>(r) 
+ \ J dr Pion (r)4>(r) - ± e leads I dr Pi (r)(f>(r) 

-kT,ie ga tes QM (6) 

where the energy levels, density, potential and induced 
charges are implicitly functions of N and the applied 
gate voltages Vi. Note that the occupation number de- 
pendence of these terms is ignored. In the T = limit 
the electrons occupy the lowest N states of the dot, and 
the free energy is denoted Fq(N, Vi). 



Z=J2 exp[-P(F({ ni }, N, V g ) - m)] 



(9) 



note that the sum on occupation configurations, {rii}, 
includes implicitly a sum on N. In Eq. |7j / is the Fermi 
function, p, is the electrochemical potential of the source 
and drain and Tp^ are the elastic couplings of level p 
to source (drain) . The notation {rii + p} denotes the set 
of occupancies {rii} with the p" 1 level, previously empty 
by assumption, filled. In Eq. M it is assumed that onl; 



single gate voltage, 
is varied. 



V g (the "plunger gate", cf. Fig. 



E. Tunneling coefficients 



The elastic couplings in Ect_ 7| are calculated from the 
self-consistent wave functions!^: 



nr np = 4 K 2 W„ 2 (a,6) 



dy f P (xb,y)x*n( x b,y) 



(10) 



where f p (xb,y) is the two dimensional part of the p th 
wave function evaluated at the midpoint of the barrier, 
Xb, and Xn( x b,y) is the n th channel wavefunction decay- 
ing into the barrier from the leads, W n (a, b) is the barrier 
penetration factor between the classical turning point in 
the lead and the point Xb, for channel n computed in 
the WKB approximation, and k is the wave vector at 
the matching point. Though the channels are ID we use 
the two dimensionaLdcnsity of states characteristic of the 
wide 2DEG regionEi 



F. Capacitance 



D. Conductance 



The master equation formula for the linear source- 
drain coad-uctance though the dot, derived by several 
authorsHaO for the case of a fixed dot spectrum, is mod- 
ified to the self-consistently determined free energy case 
as followsci: 



G{Vg) 



k B T 



p p 



Tt ' n 



x f(F({rn + P },N+ 1, V g ) - F({m}, N, V g ) - p) (7) 

where the first sum is over dot level occupation configu- 
rations and the second is over dot levels. The equilibrium 
probability distribution P eq ({rii}) is given by the Gibbs 
distribution, 

Pe q ({n t }) = ±exp[-p(F({TH},N,V g )-n)] (8) 
and the partition function is 



Quantum dot system electrostatic energies are cota- 
monly estimated on the basis of a capacitance modelEZl. 
When the self-consistent level energies and potential are 
known the total free energy can be computed without 
reference to capacitances. However, the widespread use 
of this model and the ease with which capacitances can 
be calculated from our self-consistent results (see below) 
encourages a discussion. 

For a collection of N metal elements with chargi 
and voltages Vj the capacitance matrix 

Qi = X)j=i C'u^jij can b e written i 
Green's function Gc(x,x') for Laplace's equation satis- 
fying Dirichlet boundary conditions on the element sur- 
faces: 



defined by 2 ^ 
terms of the 



C - — 

ClJ ~ 47T 2 



dSljhj ■ V x {n t ■ V a ,/G £) (x,x')) (11) 



where the integrals are over element surfaces with fij the 
outward directed normal. 
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In a system with an element of size L not much greater 
than the screening length A s , the voltage of the comj&er 
nent, and hence the capacitance, is not. well definedc3'E!3. 
In this discussed in referenced, the capacitance 

can no longer be written in terms of the solution of Pois- 
son's equation alone, but must take account of the full 
self-consistent determination of the i th charge distribu- 
tion /?i(x) from the j potential <^j(x) Vi,j. In general 
the capacitance will then become a kernel in an integral 
relation. A relationship of this kind has recently been 
derived is. terms of the Linhard screening function by 
ButtikeM 

To compute the dot self-capacitance from the calcu- 
lated self-consistent electronic structure we have three 
separate procedures. In all three cases we vary the Fermi 
energy of the dot by some small amount to change the 
net charge in the dot. This requires that the QPCs be 
closed. For the first method the total charge variation 
of the dot is divided by the change in the electrostatic 
potential minimum of the dot. This is taken as the dot 
self-capacitance Cdd- A second procedure for the dot 
self-capacitance is to divide the change in the dot charge 
simply by the fixed, imposed change of the Fermi en- 
ergy. This result is denoted C' dd . Since the change in 
the potential minimum of the dot is not always equal 
to the change of the Fermi energy these results are not 
identical. Finally, we can fit the computed free energy 
F(N, V g ) to a parabola in N at each V g . If the quadratic 
term is aN 2 then the final form for the self-capacitance 
is C' dd — l/(2a) (primes are not derivatives here). This 
form, which also serves as a consistency check on our 
functional for the energy, is generally quite close to the 
first form and we present no results for it. 

For the capacitances between dot and gates or leads, 
the extra dot charge (produced by increasing the Fermi 
energy in the dot) is screened in the gates and the leads 
so that the net charge inside the system (including that 
on the gated boundaries) remains zero. The fraction of 
the charge screened in a particular element gives that 
element's capacitance to the dot as a fraction of Cdd- 



III. RESULTS 



physical meaning. However, as pointed out by SlaterS, 
the usefulness of DF theory depends to some extent on 
being able to interpret the energies and wave functions 
as some kind of single particle spectrum. In the Coulomb 
blockade regime it is particularly important to be clear 
what that interpretation, and its limitations, are. 

A distinction is commonly made between the addi- 
tionpSnectrum and the excitation spectrum for quantum 
dotscJ'Eil. Differences between our effective single particle 
eigenvalues represent an approximation to the excitation 
spectrum. As a specific example, in the absence of de- 
polarization and excitonic effects the first single particle 
excitation from the TV-electron ground state with gate 
voltages Vi is e N+1 (N, V) - e N (N, V). 

The addition spectrum, on the other hand, depends on 
the energy difference between the ground states of the 
dot interacting with its environment at two different N. 
Thus, in our formalism, the addition spectrum is given 
by differences in F({n p }, N, Vi) at neighboring N, possi- 
bly further modulated by excitations, i.e. differences in 
the occupation numbers {n p }. 

In contrast to experiment, the electronic structure can 
be determined for arbitrary N and Vi (so long as the dot 
is closed). This includes both non-integer N as well as 
values which are far from equilibrium (differing chcxni- 
cal potential) with the leads. The "resonance curve"Q is 
given by the N which minimizes Fq(N, V g ) at each V g 
(gates other than the plunger gate are assumed fixed). 
This occurs when the chemical potential of the dot equals 
those of the leads (which are taken as equal to one an- 
other and represent the energy zero) and gives the most 
probable electron number. Results presented below as a 
function of varying gate voltage, particularly the spectra 
m Figs, landy, are assumed to be along the resonance 
curve. 



A. Electrostatics 



We consider only a small subspace of the huge available 
parameter space. For the results presented here we have 
fixed the nominal 2DEG density to 1.4 x 10 11 cm" 2 and 
the aluminum concentration of the barrier to 0.3. The 
lithographic gate pattern is shown in figure |l|, as is the 
growth profile (including our artificial second barrier). 
Some results are presented with a variation of the total 
thickness t of the AlGaAs (Fig. [l]). 

To interpret the results we note the following consider- 
ations. Hohenberg-Kohn-Sham theory provides only that 
the ground state energy of an interacting electron system 
can be written as a functional of the densityME3. The 
single particle eigenvalues e p have, strictly speaking, no 



Figure |]a shows an example of a potential profile along 
with a corresponding density plot for a quiet dot contain- 
ing 62 electrons. The basic potential/density configura- 
tion, as well as the capacitances are highly robust. These 
data are computed completely in the 2D Thomas-Fermi 
approximation, single z-subband, at T = 0.1 K. Solu- 
tion of Schrodinger's equation or variation of T result in 
only subtle changes. The depletion region spreading is 
roughly 100 nm. Figure ||b shows a set of potential and 
density profiles along the y-direction (transverse to the 
current direction) in steps of 3.3 a* B in x, from the QPC 
saddle point to the dot center. Note that the density at 
the dot center is only about 65% of the ungated 2DEG 
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FIG. 2. (a) Contour plot for density and potential, quiet 
dot, TF. Isolines in potential spaced at ~ 0.1 Ry* up to 
0.5 Ry* above Fermi level, after which much more widely. 
Density isoline spacing ~ 0.01 a* B ~ 2 , maximum density 
~ 0.1 a* B 2 - Ripples near QPCs are finite grid size effect; 
plotted x — y mesh shows every other grid line, (b) Trans- 
verse (y-direction) half-profiles of density and potential cor- 
responding to (a), taken at 3.3 a* B intervals from dot center. 
Uppermost potential trace, entirely above Fermi surface, is in 
QPC (x w 54 a* B in Fig. 2a) where density is zero. Density is 
scaled to nominal 2DEG value 0.14 a* B ~ 2 » 1.4 x 10 11 cm" 2 . 

density. Correspondingly the potential at the center is 
above the floor of the ungated 2DEG ( 0.9 Ry*). 

We discuss a simple model for the potential shape of 
a circular quantum dot below (Sec. III.B.l). Here we 
note only that the radial potential can be regarded as 
parabolic to lowest order with quartic and higher order 
corrections whose influence increases neat, the perimeter. 
In Thomas-Fermi studies on larger dotsEm with a compa- 
rable aspect ratio we find that the potential and density 



achieve only 90 % of their ungated 2DEG value nearly 
200 nm from the gate. Regarding classical, billiard cal- 
culations for gated structures therefore23~E3 even in the 
absence of impurities it is difficult to see how the "classi- 
cal" Hamiltonian at the 2DEG level can be even approx- 
imately integrable unless the lithographic gate pattern is 
azimuthally symmetricEa. 

The importance of the remote ionized impurity distri- 
bution is demonstrated in figure |^ which shows a quan- 
tum dot with randomly placed ionized 



disordered ordered, 3" = 1 /5 




FIG. 3. Contour plots of dot showing ion placement for 
disordered case (left) and ordered [T = 1/5) ion distribution, 
TF. Isolines at 0.08 Ry* up to Fermi surface, wider there- 
after. Gate voltages and locations identical in the two cases. 
Note particularly position of right QPC determined by ions 
in disordered case. 

donors on the left and with ions which have been al- 
lowed to reach quasi-equilibrium via variable range hop- 
ping, on the right. In both cases the total ion number in 
the area of the dot is fixed. The example shown here for 
the ordered case assumes, in the variable range hopping 
calculation, one ion for every five donors {T = 1/5). As 
in Ref.cl we have, for simplicity, ignored the negative U 
model for the l ip^P T | impurities (DX centers), which is still 
controversiaOO'tj. If the negative U model, at some 
barrier aluminum concentration, is correct, the most or- 
dered ion distributions will occur for T = 1/2, as opposed 
to the neutral DX picture employed hem, where ordering 
increases monotonically as T decreasesc3. 

For these assumptions figure [| indicates that ionic or- 
dering substantially reduces the potential fluctuations 
relative to the completely disordered case, even for rela- 
tively large T . Here, using ensembles of dots with varying 
T we compare the effective 2D potential with a quiet dot 
(jellium donor layer) at the same gate voltages and same 
dot electron number. The distribution of the potential 
deviation is computed as: 

P{AV) = sJpT,J2 S ( AV - ^{xi, Vj ) - V qd {x uyj )]) 

(12) 

where s labels samples (different ion distributions), typ- 
ically up to S — 10, N is the total number of x or y grid 
points in the dot (~ 50), and "qd" stands for quiet dot. 
The distributions for all T are asymmetric (Fig. ^). Al- 
though the means are indistinguishably close to zero, the 
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probability for large potential hills resulting from disor- 
der is greater than for deep depressions. Also, the distri- 
butions for points above the Fermi surface (dashed lines) 
are broader by an order of magnitude (in standard devia- 
tion) than below, due to screening. Finally, saturation as 
T -> (inset Fig. @) shows that even if the ions are ar- 
ranged in a Wigner crystal (the limiting case at T = 0) , 
potential fluctuations would be expected in comparison 
with ionic jellium. 
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FIG. 4. Histograms of deviation of effective 2D potential 
from quiet dot values at the same x,y point and same gate 
voltages, for several ion to donor ratios T, TF. Solid lines are 
statistics for points below Fermi surface, dashed lines, show- 
ing substantially more variation, above. T — 1 is completely 
random (disordered) case. Distributions uniformly asymmet- 
ric, positive potential deviations from quiet dot case being 
more likely, but means are very close to zero. Inset shows 
standard deviation of histograms versus T, triangle below, 
squares above Fermi level. 

The success of the capacitance model in describing ex- 
perimental results of charging-phenomena in mesoscopic 
systems has been remarkableEZL For our calculations as 
well, even the simplest formulations for the capacitance 
tend to produce smoothly varying results when gate volt- 
ages or dot charge are varied. Figure || shows the trend 
of the dot self-capacitances with V g . Also shown are the 
equilibrium dot electron number N and the minimum of 
the dot potential V m in as functions of V g . Note here that 
V m i n is the minimum of the 3D electrostatic potential 
rather than the effective 2D potential which is presented 
elsewhere (such as in Figs, p and ||). 

That Cdd generally decreases as the dot becomes 
smaller is-, not surprising and has been discussed 
elsewhereEJ. All three forms of Cdd are roughly in agree- 
ment giving a value ~ 2 fF (the capacitance as calcu- 
lated from the free energy is not shown) . The fluctuations 
result from variations in the quantized level energies as 
the dot size and shape are changed by V g . Note that 
numerical error is indiscernible on the scale of the fig- 



ure. The pronounced collapse of C' dd near V g — —1.15 V, 
which is expanded in the upper panel, shows the presence 
of a region where the change of N with Ep is greatly 
suppressed. Since the change of V m in with Ep is sim- 
ilarly suppressed there is no corresponding anomaly in 
Cdd- Interestingly, the capacitance computed from the 
free energy also reveals no deep anomaly. 
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FIG. 5. Dot self-capacitances, equilibrium electron number 
and potential minimum as a function of plunger gate voltage 
(lower). Numerical uncertainty is indiscernible, so variations 
of Cdd are real and related to spectrum. C dd calculated using 
AEf rather than AVmin, so strong anomaly near —1.15 V 
due to rigidity of N. Upper panel: expanded view of capaci- 
tances near anomaly; cf. spectrum, Fig. 9. 

The anomaly at V g — —1.15 V and also the fluctuation 
in the electrostatic properties near —1.1 V are related to 
a shell structure in the spectrum which we discuss below. 

A frequently encountered model for the classical charge 
distribution in a quantum dot is the circular_cpnduct- 
ing disk with a parabolic confining potentialli-lcj. It can 
be shown (solving, for example, Poisson's equation in 
oblate spheroidal coordinates) that for such a model the 
2D charge distribution in the dot goes as 



n{r) = ? i(0)(l-r 2 /i? 2 ) 1/2 



(13) 



where R is the dot radius and n(0) = 3N/2irR 2 is the 
density at the dot center. The "external" confining po- 



tential is assumed to go as V(r) = Vo 
related to N through 



R 



^£1 

4 nk 



N 



kr 2 /2 and R is 



(14) 



where k is the dielectric constantB. _ 

To justify this model, the authors of Ref.Q claim that 
the calculations of Kumar et alx3 show that "the con- 
finement... has a nearly parabolic form for the external 
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confining potential (sic)." This is incorrect. What Ku- 
mar et al.'s calculations shows is that (for N ~ 12) 
the self-consistent potential, which includes the poten- 
tial from the electrons themselves, is approximately cut- 
off parabolic—-. The external confining potential, as it is 
used in Refic-a, would be that produced by the donor 
layer charge and the charge on the surface gates only. 
We introduce a simple model (see III.B.l below) wherein 
this confining potential charge is replaced by a circular 
disk of positive charge whose density is fixed by the dop- 
ing density and whose radius is determined by the num- 
ber of electrons in the dot. The gates can be thought 
of as merely cancelling the donor charge outside that ra- 
dius. The essential point, then, is this: adding electrons 
to the dot decreases the (negative) charge on the gates 
and therefore increases— the radius. One can make the 
assumption, as in Ref.E-3, that the external potential is 
parabolic, but it is a mistake to treat that parabolicity, 
k, as independent of N. 

This is illustrated in figure ^ where we have plotted 
contours for the change in the 2D density, as Ep is in- 
crementally increased, 
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FIG. 6. Grey scale of density change as Fermi energy in dot 
is raised relative to leads, Thomas- Fermi (TF). Total change 
in TV about 1.4 electrons. Screening charge, white region, in 
leads is positive. White curve gives profile along line bisecting 
dot, scaled to average change of TV" per unit area. Right panel 
shows model of Ref. 44 where confining potential has fixed 
parabolicity. Note that this model drastically underestimates 
degree to which charge is added to perimeter. 

as determined self-consistently (Thomas-Fermi every- 
where, left panel) and as determined from Eq. [l3| The 
white curves display the density change profiles across 
the central axis of the dot. The total change in N is 
the same in both cases, but clearly the model of Eq. [l3| 
underestimates the degree to which new charge is added 
mostly to the perimeter. 

Recently the question of charging energy renormaliza- 
tion via tunneling as the conductance Go through-a QPC 
approaches unity has received much attentiorO" ufj. In a 
recent experiment employing two dots in series a splitting 
of the Coulomb oscillation peaks has been observed as the 
central QPC (between the two dots) is loweredEJ. Per- 
turbation theory for small Go and a model which treats 



the decaying channel between the dots as a Luttinger liq- 
uid for Go — ► 1 (e 2 /h) lead to expressions for the peak 
splitting which is linear in Go in the former case and goes 
as (1 — Go)ln(l — Go) in the latter case. 
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FIG. 7. Variation of dot capacitances with QPC voltage 
(gates 1 and 4 in figure 1). Solid lines for Vum are ef- 
fective 2D potential for left (right) saddle point (right hand 
scale). Cs{A) and Cs(B) are dot self-capacitances (cf. Fig. 
5) computed using AVmin and AEf respectively. "Source" 
is (arbitrarily) outside left saddle point. Note that Vl goes 
practically to zero but the dot capacitance to the source only 
marginally increases relative to dot to drain capacitance. Ca- 
pacitances for QPC and plunger are for a single finger only in 
each case. Anomaly related to dot reconstruction also visible 
here as QPC voltage is changed. 

A crucial assumption of the model, however, is that 
the "bare" capacitance, specifically that between the dots 
Cdi-d2, remains approximately independent of the height 
of the QPC, even when an open channel connects the 
two dots. Thus the mechanism of the peak splitting is 
assumed to be qualitatively different from a model which 
predicts peak splitting entirely on an electrostatic basis 
when the inter-dot capacitance increases greatlyc2l. The 
independence of Cdi-d2 from the QPC potential is plau- 
sible insofar as most electrons, even when a channel is 
open, are below the QPC saddle points and hence lo- 
calized on either one dot or the other. Further, if the 
screening length is short and if the channel itself does 
not accommodate a significant fraction of the electrons, 
there is little ambiguity in retaining Cd\-d2 to describe 
the gross electrostatic interaction of the dots, even when 
the dots are connected at the Fermi level. 

In figure we present evidence for this theory by show- 
ing the capacitance between a dot and the leads as the 
QPC voltage is reduced. In the figure Vum is the effec- 
tive 2D potential of the left (right) saddle point as the 
left QPC gate voltages Vqpc only are varied. The dot 
is nearly open when the QPC voltages (both pins on the 
left) reach ~ —1.34 V. The results here use the full quan- 
tum mechanical solution (without the LDA exchange- 
correlation energy), however the electrons in the lead 







continue to be treated with a 2D TF approximation. The 
dot "reconstruction" seen in figure H is visible 
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FIG. 8. AlGaAs thickness dependence of capacitances 
(lower). Self-capacitance decreases as gates get closer to 
2DEG. Upper panel shows that, for smaller t the poten- 
tial confinement is steeper and charge more compact, hence 
smaller Cdd- t\ = 5.25, t% = 7.5, t$ = 9, and = 12a* B . Rel- 
ative capacitance from dot to gates and leads fairly insensitive 
to t. 

here also around Vqpc = —1.365 V. Note that the 
right saddle point is sympathetically affected when we 
change this left QPC. While the effect is faint, ~ 5% 
of the change of the left saddle, the sensitivity of tun- 
neling to saddle point voltage (see also below) has re- 
sulted in this kind of cross-talk being problematical for 
experimentalists. The figure also shows that the capac- 
itance between the dot and one lead exceeds that to a 
(single) QPC gate or even to a plunger gate. However 
the most important result of the figure is to show that 
the dot to lead capacitance is largely insensitive to QPC 
voltage. When the left QPC is as closed as the right 
(Vqpc ~ —1.375 V) the capacitances to the source and 
drain are equal. But even near the open condition the ca- 
pacitance to the left lead (arbitrarily the "source") only 
exceeds that to the drain (which is still closed) minutely. 
Therefore the assumptions of a "bare" capacitance which 
remains constant even as contact is made with a lead (or, 
in the experiment, another dot) seems to be very well 
founded. 

As noted above, the interaction between a gate and the 
2DEG depends upon the distance of the gates from the 
2DEG, i.e., the AlGaAs thickness t. In figure ^ we show 
that, as we decrease t, simultaneously changing the gate 



voltages such that N and the saddle point potentials re- 
main constant, the total dot capacitance also decreases, 
but the distribution of the dot capacitance between leads, 
gates and (not shown) back gate change only moderately. 
That gates closer to the 2DEG plane should produce dots 
of lower capacitance is made clear in the upper panel of 
the figure, which shows the potential and density profile 
(using TF) near a depletion region at the side of the dot 
at varying t and constant gate voltage. For smaller t the 
depletion region is widened but the density achieves its 
ungated 2DEG value (here 0.14 a* B ~ 2 ) more quickly; a 
potential closer to hard walled is realized. In the pres- 
ence of stronger confinement the capacitance decreases 
and the charging energy increases. 

The profile of the tunnel barriers and the barrier pen- 
etration factors are also dependent on t. However we 
postpone a discussion of this until the section on tunnel- 
ing coefficients. 



B. Spectrum 

The bulk electrostatic properties of a dot are, to first 
approximation, independent of whether a Thomas-Fermi 
approximation is used or Schrodinger's equation is solved. 
A notable exception to this is the fluctuation in the ca- 
pacitances. Figure |j] shows the plunger gate voltage de- 
pendence of the energy levels. The Fermi level of the 
dot is kept constant and equal to that of the leads (it 
is the energy zero). Hence as the gate voltage increases 
(becomes less negative) N increases. 

Since the QPCs lie along the x-axis, the dot is never 
fully symmetric with respect to interchange of x and 
y, however the most symmetric configuration occurs for 
V g ~ —1.16 V, towards the right side of the plot. The 
levels clearly group into quasi-shells with gaps between. 
The number of states per shell follows the degeneracy of 
a 2D parabolic potential, i.e. 1,2,3,4,... degenerate levels 
per shell (ignoring spin) . There is a pronounced tendency 
for the levels to cluster at the Fermi surface, here given 
by E = 0, which we discuss below. 



1. Shell structure 

Shell structure in atoms arises from the approximate 
constancy of individual electron angular momenta, and 
degeneracy with respect to z-projection. Since in two 
dimensions the angular momentum m is fixed in the z 
(transverse) direction, the isotropy of space is broken and 
the only remaining manifest degeneracy, and this only for 
azimuthally symmetric dots, is with respect to ±z. A two 
dimensional parabolic potential, in the absence of mag- 
netic field, possesses an accidental degeneracy for which 
a shell structure is recovered. 

We have shown above that modelling a quantum dot 
as a classical, conducting layer in an external parabolic 
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potential kr 2 /2, where k is independent of the number of 
electrons in the dot, ignores the image charge in the sur- 
face gates forming the dot and therefore fails to properly 
describe the evolving charge distribution as electrons are 
added to the dot. A more realistic model, which explains 
the approximate parabolicity of the self-consistent poten- 
tial, and hence the apparent shell structure, is illustrated 
in figure [l^. 



A simple calculation for the radial potential (for a < 
R) in the electron layer (z = 0) gives, for the first few 
terms: 



_ 15 a 4 r 2 i 45_ a 2 r 4 , 
32 WW ' 128 WW " t " 



— - 

sWW 



(15) 



disorder order 




I 1 1 1 1 I 
-1.1 -1.15 
Gate Voltage V, [VI 

FIG. 9. Electronic spectrum showing level grouping into 
shells for quiet dot (Hartree), quiet dot with LDA ex- 
change-correlation, disordered sample !F = 1 and ordered 
sample T = 1/5. Range of gate voltage in latter three is 



from V„ 



-1.142 to -1.17 V. 



The basic electrostatic structure of a quantum dot, in 
the simplest approximation, can be represented by two 
circular disks, of radius R and homogeneous charge den- 
sity (To, separated by a distance a. The positive charge 
outside R is assumed to be cancelled by the surface gates. 
This approximation will be best for surface gates very 
close to the donor layer (i.e. small t). Larger AlGaAs 
thicknesses will require a non-abrupt termination 




FIG. 10. Schematic for simple two charge disk model of 
quantum dot. Positive charge outside radius R taken to be 
uniformly cancelled by gates, electric charge in 2DEG mirrors 
positive charge. Resultant radial potential in 2DEG plane, 
Eq. 15, dominated by parabolic term inside R. 

of the positive charge. In cither case, the electronic 
charge is assumed in the classical limit to screen the back- 
ground charge as nearly as possible. This is similar to the 
postulate in which wide parabolic quantum wells are ex- 
pected to producpiapproximately homogeneous layers of 
electronic chargeEi 



where Ne = ■kR 2 oq and k is the background dielectric 
constant. While the coefficient of the quartic term is com- 
parable to that of the parabolic term, the dependences 
are scaled by the dot radius R. Hence, the accidental 
degeneracy of the parabolic potential is broken only by 
coupling via the quartic term near the dot perimeter. 
This picture clearly agrees with the full self-consistent 
results wherein the parabolic degeneracy is observed for 
low lying states and a spreading of the previously degen- 
erate states occurs nearer to the Fermi surface. 

Comparison (not shown) of the potential computed 
from Eq. 
Fig. |b) 



151 and the radial potential profile (lowest curve, 
from the full self-consistent structure, shows 
good agreement for overall shape. However the former is 
about 25% smaller (same N) indicating that the sharp 
cutoff of the positive charge is, for these parameters, too 
extreme. However Eq. |l5| improves for larger N and/or 
smaller t. 

The wavefunction moduli squared associated with the 

Fig. g quiet dot levels for V g 1.16 V, N w 54 are 

shown schematically for levels 1 through 10 in figure O, 
and for levels 11 through 35 in figure [jjl 

The lowest level in a shell is, for the higher shells, typi- 
cally the most circularly symmetric. When the last mem- 
ber of a shell depopulates with V g the inner shells expand 
outward, as can be seen near V g = —1.15 V (Fig. g) 
where level p = 29 depopulates. Since to begin filling a 
new shell requires the inward compression of the other 
shells and hence more energy, the capacitance decreases 
in a step when a shell is depopulated. The shell structure 
should have two distinct signatures in the standard (elec- 
trostatic) Coulomb oscillation experimentE3. First, since 
the self-capacitance drops appreciably (figure^) when the 
last member of a shell depopulates, here N goes from 57 
to 56, a concomitant discrete rise in the activation energy 
in the minimum between Coulomb oscillations can b|& 
predicted. Second, envelope modulation of peak heightscl 
occurs when excited dot states are thermally accessible 
as channels for transport, as opposed to the T — case 
where the only channel is through the first open state 
above the Fermi surface (i.e. the N + 1 st state). When 
N is in the middle of a shell of closely spaced, spin de- 
generate levels, the entropy of the dot, fcsZnfi, where f2 
is the number of states accessible to the dot, is sharply 
peaked. For example, for six electrons occupying six spin 
degenerate levels (i.e. twelve altogether) 
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FIG. 11. Schematic showing the first ten levels of quiet dot. 
Shell structure consistent with n + m = constant, where n and 
m are nodes in x and y. Lower energy states show rectangular 
symmetry. 
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FIG. 12. Levels 11 through 35 (each spin degenerate) of 
quiet dot, Hartree. Circular symmetry increases with energy. 
States elongated in x (horizontal) most connected to leads. 

all within ksT of the Fermi surface, the number of 
channels available for transport is 924. For eleven elec- 
trons in the shell, however, the number of channels re- 
duces to 12. Consequently, minima and peaks of envelope 
modulation (see also figure ^l] below) of CB oscillations 
which are frequently observed are clear evidence of level 
bunching, if not an organized shell structure. 

Recently experimental evidence has accumulated for 
the existence of-a shell structure as observed by inelastic 
light scattering^ and via Coulomb oscillation peak posi- 
tions in transport through extremely small (N ~ — 30) 
vertical quantum dotsEj. Interestingly, a classical treat-, 
ment, via Monte-Carlo molecular dynamics simulationE-3 
also predicts a shell structure. Here, the effect of the neu- 
tralizing positive background are assumed to produce a 
parabolic confining potential. A similar assumption is 
made in Ref.E3 which analyzes a vertical structure simi- 
lar to that of Reffi We believe that continued advances 
in fabrication will result in further emphasis on such in- 
variant, as opposed to merely statistical, properties of 
dot spectra. 

As noted above, there is a strong tendency for levels 
at the Fermi surface to "lock." Such an effect has been 



described by Sun et alxB in the case of subband levels for 
parallel quantum wires. In dots, the effect can be viewed 
as electrostatic pressure on the individual wavefunctions 
thereby shifting level energies in such a way as to produce 
level occupancies which minimize the total energy. Inso- 
far as a given set of level occupancies is electrostatically 
most favorable, level locking is a temperature dependent 
effect which increases as T is lowered. This self-consistent 
modification of the level energies can also be viewed as 
an excitonic correction to excitation energies. 

The difference between the cases of a quantum dot 
and that of parallel wires is one of localized versus ex- 
tended systems. It is well known that, unlike Hartree- 
Fock theory, wherein self-interaction is completely can- 
celled since the direct and exchange terms have the same 
kernel l/|r — r'|, in Hartree theory and even density func- 
tional theory in the LDA, uncorrected self-interaction 
remainsE3. While it is reasonable to expect that ex- 
cited states will have their energies corrected downward 
by the remnants of an excitonic effect, we expect that 
LDA and especially Hartree calculations will generally 
overestimate this tendency to the extent that corrections 
for self-interaction are not complete. 

The panel labelled "xc" in figure || illustrates the pre- 
ceding point. In contrast to the large panel (on the left) 
these results have had the XC potential in LDA included. 
The differences between Hartree and LDA are generally 
subtle, but here the clustering of the levels at the Fermi 
surface is clearly mitigated by the inclusion of XC. The 
approximate parabolic degeneracy is evidently not bro- 
ken by LDA, however, and the shell structure remains in- 
tact. Similarly for xc, the capacitances also show anoma- 
lies near the same gate voltages, where shells depopulate, 
as in figure S, which is pure Hartree. 

The two remaining panels in figure || illustrate the ef- 
fects of disorder and ordering in the donor layer (XC not 
included). As with the "xc" panel, V g is varied between 
— 1.142 and —1.17 V. The "disorder" panel represents 
a single fixed distribution of ions placed at random in 
the donor layer as discussed above. Similarly, the "or- 
der" panel represents a single ordered distribution gen- 
erated froBif-a random distribution via the Monte-Carlo 
simulationO'EJ; here T = 1/5 (cf. two panels of Fig. ||). 

The shell structure, which is completely destroyed for 
fully random donor placement (see also Fig. 14), is al- 
most perfectly recovered in the ordered case. In both 
cases the energies are uniformly shifted upwards relative 
to the quiet dot by virtue of the discreteness of donor 
charge (cf. also discussion of Fig. [| above). Closer ex- 
amination of the disordered spectrum shows considerably 
more level repulsion than the other cases. 

The application of a small magnetic field, roughly a 
single flux quantum through the dot, has a dramatic im- 
pact on both the spectrum, figure |l3|, and the wave func- 
tions, figure |lj, top. The magnetic field dependence of 
the levels (not shown) up to 0.1 T exhibits shell splitting 
according to azimuthal quantum number as 
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FIG. 13. V g dependence at fixed B (0.05 T) of level ener- 
gies, quiet dot. Multiple re-constructions seen as levels de- 
populate. Homogeneous level spacing related to uniformity 
of Coulomb oscillation peak heights in a magnetic field. 

well as level anti-crossing. By 0.05 T level spacing 
(Fig. 13) is substantially more uniform than B = 0, 
Fig. |. Furthermore, while the B = quiet dot dis- 
plays reconstruction due to the depopulation of shells at 



V„ 



-1.15 and -1.1 V, the B = 0.05 T results show a 



similar pattern, a step in the levels, repeated 
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FIG. 14. Levels 31 through 35 for (from bottom) quiet dot 
with LDA for XC, Hartree for disordered dot, Hartree for or- 
dered dot T = 1/5 and B = 0.05 T. XC changes ordering of 
some levels, but has very little influence on states. Ordered 
case recovers much of quiet dot symmetry. Small B changes 
states altogether. 

many times in the same gate voltage range. The phys- 
ical meaning of this is clear. The magnetic field princi- 



pally serves to remove the azimuthal dependence of the 
mod squared of the wave functions (Fig. |TJ). In a mag- 
netic field, the states at the Fermi surface also tend to 
be at the dot perimeter. Depopulation of an electron in 
a magnetic field, like depopulation of the last member 
of a shell for B = 0, therefore removes charge from the 
perimeter of the dot and a self-consistent expansion of 
the remaining states outward occurs. 



C. Statistical properties 

1. Level spacings 

The statistical spectral properties of quantum systems 
whose classical Hamiltonian is chaotic are believed rta 
obey the predictions of random matrix theory (RMT)l3. 
Arguments for this conjecture however invariably treat 
the Hamiltonian as a large finite matrix with averaging 
taken only near the band center. Additionally, an often 
un-clearly stated assumption is that the system in ques- 
tion can be treated semi- classically, that is, in some sense 
the action is large on the scale of Planck's constant and 
the wavelength of all relevant states is short on the scale 
of the system size. Clearly, for small quantum dots these 
assumptions are violated. 

RMT predictions apply to level spacings S and to 
transition -amplitudes (for the "exterior problem," level 
widths r)Ej. RMT is also applied to scattering ma- 
trices in ixtvestigations of transport properties of quan- 
tum wirea£j. Ergodicity for chaotic systems is the claim 
that variation of some external parameter X will sweep 
the Hamiltonian rapidly through its entire Hilbert space, 
whereupon energy averaging and ensemble (i.e. X) av- 
eraging produce identical statistics. In our study X is 
either the set of gate voltages, the magnetic field or the 
impurity configuration and we consider the statistics of 
the lowest lying 45 levels (spin is ignored here). Care 
must also be taken in removing the secular variations of 
the spacings or widths with energy, the so-called unfold- 
ing. 

According to RMT level repulsion leads to statistics of 
level spacings which are given by the "Rayleigh distribu- 
tion:" 



P(S) 



irS 
2D 



exp(-nS 2 /4D 2 



(16) 



where D is the mean local spacin ME. Figure ^| shows 
the calculated histogram for the level spacings for the 
quiet dot as well as for disordered, ordered and ordered 
with B — 0.05 T cases. Statistics are generated from 
(symmetrical) plunger gate variation, in steps of 0.001 V, 
over a range of 0.1 V, employing the spacings between the 
lowest 45 levels; thus about 4500 data points. Deviation 
from the Rayleigh distribution is evident. An important 
feature of our dot is symmetry under inversion through 
both axes bisecting the dot. It is well known that groups 
of states which are 



13 



1.0 



0.5 



0.0- 
0.5- 



\ ordered 

iV 57=1/5 



disordered 




v B=0.05 T 



1 2 3 4 5 

level spacing S/D 

FIG. 15. Histograms of level spacings, normalized to local 
level spacing. Dark curve represents Rayleigh distribution. 
Black bars (main panel) include all states, white bars only 
for states that are completely even under x or y inversion. 
Insets: disordered panel recapitulates Rayleigh distribution, 
both ordered and B ^ marginally but significantly different. 



un-coupled will, when plotted together, show a Poisson 
distribution for the spacings rather than the level repul- 
sion of Eq. [l(]. Thus we have also plotted (white bars) 
the statistics for those states which are totally even in 
parity. While the probability of degeneracy decreases, a 
X 2 test shows that the distribution remains substantially 
removed from the Rayleigh form. 

In contrast to this, the disordered case shows remark- 
able agreement with the RMT prediction. As with the 
spectrum in figure ^ we use a single ion distribution. 
However we also find (not shown) that fixing the gate 
voltage and varying the random ion distributions results 
in nearly the same statistics. When the ions are allowed 
to order the level statistics again deviate from the RMT 
model. This is somewhat surprising since Fig. || shows 
that, even for T = 1/5, the standard deviation of the 
effective 2D potential below the Fermi surface from the 
quiet dot case, ~ 0.05 Ry*, is still substantially greater 
than the mean level spacing ~ 0.02 Ry* . We have re- 
cently shown that, as T goes from unity to zero, a con- 
tinuous transition from the level repulsion ofiiq. |l^ to a 
Poisson distribution of level spacings resultsEx Finally, 
the application of a magnetic field strong enough to break 
time-reversal symmetry clearly reduces the incidence of 
very small spacings, but the distribution is still signifi- 
cantly different from RMT. 



the barrier, P n = W n (a,c) where c is the classical turn- 
ing point on the dot side of the barrier, is plotted as a 
function of QPC voltage in figure [l6|. P n is simply the 
WKB penetration for a given channel with a given self- 
consistent barrier profile, and can be computed at any 
energy. Here we have computed it at energies coincident 
with the dot levels. Therefore the dashes recapitulate the 
level structure, spaced now not in energy but in "bare" 
partial width. The actual width of a level depends upon 
the wave function for that state (cf. Eq. [l(]). For ener- 
gies above the barrier ln(P) = 0. The solid lines repre- 
sent P at the Fermi surface computed for three different 
AlGaAs thicknesses t (as in figure ||) and for both n = 1 
and n — 2 (the dashes are computed for t = 12 a* B ). The 
QPC voltage is given relative to values at which P for 
n = 1 is the same for all three t (hence the top three 
solid lines converge at AVqpc = 0). 

Quite surprisingly t has very little influence on the 
trend of P with QPC voltage. Note that the ratio of 
barrier penetration between the second and first chan- 
nels P2/P1 decreases substantially with increasing t since 
the saddle profile becomes wider for more distant gates. 
Even for t = 7.5 a* B however, penetration via the second 
channel is about a factor of five smaller than via n — 1. 
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FIG. 16. Barrier penetration factors from classical turning 
point in lead to turning point in dot at same energy, as a 
function of QPC voltage offset. P evaluated at energies of 
states in quiet dot for AlGaAs thickness t = 12 a* B . Solid 
lines indicate barrier penetration at Fermi level. Upper three 
lines for first channel, t — 7.5, 9.0, 12.0 a* B respectively. Lower 
three lines for second channel, same t. AVqpc zero set such 
that first channel conducts equally at the Fermi surface for 
all t. 



2. Level widths 

In Eq. [l(] we defined W n (a, b) as the barrier pene- 
tration factor from the classically accessible region of 
the lead to the matching point in the barrier, for the 
n th channel. The penetration factor completely through 



Figure |l7| shows the partial width for tunneling via 
n = 1 through the barrier, now using the full Eq. |l0|, for 
the quiet dot. The barriers here are fairly wide. While 
this strikingly coherent structure is quickly destroyed by 
discretely localized donors even when donor ordering is 
allowed, the pattern is nonetheless highly informative. 
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The principal division between upper and lower states is 
based on parity. States which are odd with respect to 
the axis bisecting the QPC should in fact have identi- 
cally zero partial width (that they don't is evidence of 
numerical error, mostly imperfect convergence). 
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FIG. 17. Partial widths (through first channel) for tunnel- 
ing to the leads, quiet dot. Numbers indicate ordinate of wave 
functions, Figs. 11 and 12. Weakly connected states zero by 
parity (non-zero only through numerical error). 

Note that this division is largely preserved for discrete 
but ordered ions. The widest states (largest T) are la- 
belled with their level index for comparison with their 
wave functions in Figs. [l2| and [13]. Comparison shows 
they represent the states which are aligned along the di- 
rection of current flow. Thus in each shell there arc likely 
to be a spread of tunneling coefficients, that is, two mem- 
bers of the same shell will not have the same T. 

Statistics of the level partial widths are shown in figure 
|l8| , here normalized to their local mean values. While the 
statistics for the quiet dot are in substantial disagreement 
with RMT it is clear that discreteness of the ion charge, 
even ordered, largely restores ergodicity. The RMT pre- 
diction, the "Porter-Thomas" (PT) distribution, is also 
plotted. For non-zero B, panels (b) and (c), the predicted 
distribution is x! rather than PT. Even the completely 
disordered case (e) retains a fraction of vanishing par- 
tial width states. Since in our case the zero width states 
result from residual reflection symmetry, it would be in.- 
teresting to compare the data from references^] ando, 
which employ nominally symmetric and non-symmetric 
dots respectively, to see if the incidence of zero width 
states shows a statistically significant difference. 

One further statistical feature which we calculate is the 
autocorrelation function of the level widths as an external 
parameter X is varied: 



C(AX) 



STi(Xj)STi(Xj + AX) 



Eij 5r i (x J rJ^~sr~{x~TAxy^ 



(17) 



where ST^X) = T l (X)-T l (X), and where T(X) is again 
the local average, over levels at fixed X, of the level 
widths. Note that the sum on i is over levels and the 
sum on j is over starting values of X. 




FIG. 18. Statistics of unfolded partial level widths, first 
channel only, (a) quiet dot showing large weight near zero 
due to parity, (b) and (c) have B — 0.05 T, quiet dot and 
disordered, respectively. Remnant of peak at small coupling 
remains. Dark line represents \2 distribution predicted by 
RMT. (d) and (e) are ordered and disordered with B = 0. 
Ordered case differs significantly from Porter-Thomas distri- 
bution plotted in black here. 

In figure [l^ we show the autocorrelation function for 
varying magnetic field (cf. Refo, figure |||). The sample 
is ordered, T = 1/5. 
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FIG. 19. Autocorrelation function for level partial widths; 
ordered, T = 1/5, averaged over B starting point and all 45 
levels. Range of B is only — 0.1 T, so statistics are weaker 
to the right. Pronounced anti-correlation near 0.03 T in con- 
tradiction with RMT. 
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Our range of B only encompasses [0,0.1] T in steps 
of 0.005 T, so we have here averaged over all levels (i.e. 
i = 1 — i-45). The crucial feature, which has been noted 
inJjlefsO and, for conductance correlation in open dots 
ir£3, is that the correlation function becomes negativejn 
contradiction with a recent prediction based on RMTo. 
Indeed, as noted by Bird et alx3, an oscillatory structure 
seems to emerge in the data. Comparison with calcula- 
tion here is hampered since the statistics are less good as 
B increases. 

Nonetheless, the RMT prediction is clearly erroneous. 
We speculate th^t the basis of the discrepancy is in 
the assumntk.r£3 that C(AX) = C(-AX). Given this 
assumptiorEJ the correlation becomes positive definite. 
Physically this means that, regardless of whether B is 
positive or negative, the self-correlation of a level width 
will be independent of whether AB is positive or nega- 
tive. This implies that the level widths should be inde- 
pendent of the absolute value of B, or any even powers of 
B, at least to lowest order in AB/B. For real quantum 
dot systems this assumption is inapplicable. 

Similar behaviour is observed with X taken as the 
(plunger) gate voltage, for which we have considerably 
more calculated results, Fig. EG. 




AVJV) 

FIG. 20. Autocorrelation function with V g , averaged over 
groups of 15 levels (upper panel). Number indicates center of 
(contiguous) range of averaged values. Dashed line is average 
of all states. Lower panel is grey scale for autocorrelation of 
individual levels averaged only over V g starting point. Black 
is 1.0 and white is —1.0. Data suggests that behaviour of 
autocorrelation is sensitive to which levels are averaged. 

The upper panel is the analogue of Fig. [H], only we 
have broken the average on levels into separate groups 
of fifteen levels centered on the level listed on the figure 
(e.g., the "28" denotes a sum in equ. 0of i = 21,35). 
the lower panel shows the autocorrelation as a grey scale 
for the individual levels (averaging performed only over 
starting V g ). The very low lying levels, up to ~ 10, re- 



main self-correlated across the entire range of gate volt- 
age. This simply indicates that the correlation field is 
level dependent. However, rather than becoming uni- 
formly-, grey in a Lorentzian fashion, as predicted by 
RMTE3, individual levels tend to be strongly correlated 
or anti-correlated with their original values, and the dis- 
appearance of correlation only occurs as an average over 
levels. 

Again we expect that the explanation for this be- 
haviour lies in the shell structure. Coulomb interaction 
prevents states which are nearby in energy from having 
common spatial distributions. Thus in a given range of 
energy, when one state is strongly connected to the leads, 
other states are less likely to be. Further, the ordering 
of states appears to survive at least a small amount of 
disorder in the ion configuration. 



D. Conductance 



The final topic we consider here is the Coulomb oscil- 
lation conductance of the dot. We will here focus on the 
temperature dependences, although statistical properties 
related to ion ordering are also interesting. 

We have shown in Ref.l that detailed temperature de- 
pendence of Coulomb oscillation amplitudes can be em- 
ployed as a form of quantum dot spectroscopy. Roughly, 
in the low T limit the peak heights give the individual 
level connection coefficients and, as temperature is raised 
activated conductance at the peaks depends on the near- 
est level spacings at the Fermi surface. In this regard 
we have explained envelope modulation of peak heights, 
which had previously not been understood, as clear evi- 
dence of thermal activation involving tunneling through 
excited states of the dota. 

Figure |2l| a shows the conductance as a function of 
plunger gate voltage for the ordered dot at T = 250 mK. 
Note that the magnitude of the conductance is small be- 
cause the coupling coefficients are evaluated with rela- 
tively wide barriers for numerical reasons. Over this 
range the dot N depopulates from 62 (far left) to 39. The 
level spacings and tunneling coefficients are all changing 
with V g . At low temperature a given peak height is deter- 
mined mostly by the coupling to the first empty dot level 
(Tjv+i) an d by the spacings between the N th level and 
the nearest other level (above or below) . The relative im- 
portance of the T's and the level spacings can obviously 
vary. In this example, Figs. |2l| a and [2l| b suggest that 
peak heights correlate more strongly with the level spac- 
ings. The double envelope coincides with the Fermi level 
passing through two shells. In general, the DOS fluctua- 
tions embodied in the shell structure and the observation 
(above) that within a shell a spreading of the T's (with 
a most strongly coupled level) results from Coulomb in- 
teraction provide the two fundamental bases of envelope 
modulation. 
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FIG. 21. (a) Conductance versus V g for ordered dot, 
T — 0.25 K. (b) Fermi surface level spacing and tunnel- 
ing coefficient at resonance. Conductance in (a) correlates 
somewhat more strongly with smaller level spacing than with 
larger V. 

Finally, we typically find that, when peak heights are 
plotted as a function of temperature (not shown) some 
peaks retain activated conductance down to T = 10 mK. 
Since the dot which we are modelling is small on the 
scale of currently fabricated structures, this study sug- 
gests that claims to have reached the regime where all 
Coulomb oscillations represent tunneling through a sin- 
gle dot level are questionable. 



IV. CONCLUSIONS 



We have presented extensive data from calculations 
on the electronic structure of lateral GaAs — AlGaAs 
quantum dots, with electron number in the range of 
N = 50 — 100. Among the principal conclusions which 



we reach are the following. 

The electrostatic profile of the dot is determined by 
metal gates at fixed voltage rather than a fixed space 
charge. As a consequence of this the model of the dot as 
a conducting disk with fixed, "external," parabolic con- 
finement is incorrect. Charge added to the dot resides 
much more at the dot perimeter than this model pre- 
dicts. 

The assumption of complete disorder in the donor layer 
is probably overly pessimistic. In such a case the 2DEG 
electrostatic profile is completely dominated by the ions 
and it is difficult to see how workable structures could be 
fabricated at all. The presence of even a small degree of 
ordering in the donor layer, which can be experimentally 
modified by a back gate, dramatically reduces potential 
fluctuations at the 2DEG level. 

Dot energy levels show a shell structure which is robust 
to ordered donor layer ions, though for complete disorder 
it appears to break up. The shell structure is responsi- 
ble for variations in the capacitance with gate voltage 
as well as envelope modulation of Coulomb oscillation 
peaks. The claims that Coulomb oscillation data through 
currently fabricated lateral quantum dots shows unam- 
biguous transport through single levels are questionable, 
though some oscillations will saturate at a higher tem- 
perature than others. 

The capacitance between the dot and a lead increases 
only very slightly as the QPC barrier is reduced. Thus 
the electrostatic energy between dot and leads is dom- 
inated by charge below the Fermi surface and splitting 
of oscillation peaks through double dot structured^ is 
undoubtedly a result of tunneling. 

Finally, chaos is well known to be mitigated in 
quantum .-systems where barrier penetration is non- 
negligibleE3. Insofar as non-inegrability of the underlying 
classical Hamiltonian is being-used as the justification for 
an assumption of ergodicityEj in quantum dots, our re- 
sults suggest that further success in comparison with real 
(i.e. experimental) systems will occur only wliea,account 
is taken in, for example, the level velocityE3iZ3, of the 
correlating influences of quantum mechanics. 
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